/*
 * td2.cpp
 *
 *  Created on: 29 janv. 2012
 *      Author: guillaume Bouhier, anne-elise Poirot
 */


#include "td2.hpp"

void analyse_97_lifting(double *x, int p, FILE* file){

	assert((((p & -p) xor p) == 0)  &&  (p != 0) && "P ,n'est pas une puissance de 2");
	double a = 0;
	double* temp = new double[p];
	int deux_n_plus_un, deux_n_plus_deux, deux_n_moins_un;
	for(int n=0; n < p/2; n++){

		//deux_n_plus_un =  2*n+1;
		//if(deux_n_plus_un > p-1) ind -= p;

		//Pas d'effet de bord dans l'AMR2D avec cette methode de bouclage
		deux_n_plus_deux = 2*n+2;
		if(deux_n_plus_deux > p-1) deux_n_plus_deux -= p;

		deux_n_moins_un =  2*n-1;
		if(deux_n_moins_un < 0) deux_n_moins_un += p;

		//Methode mirroir
		/*deux_n_plus_deux = 2*n+2;
		if(deux_n_plus_deux > p-1) deux_n_plus_deux -= (deux_n_plus_deux-p-1);

		deux_n_moins_un =  2*n-1;
		if(deux_n_moins_un < 0) deux_n_moins_un = -deux_n_moins_un;*/


		//Prediction 1
		a = -1.586134342;
		x[2*n+1] = a*x[2*n] + x[2*n+1] + a*x[deux_n_plus_deux];
	}

	for(int n=0; n < p/2; n++){

		deux_n_moins_un =  2*n-1;
		if(deux_n_moins_un < 0) deux_n_moins_un += p;

		//MaJ 1
		a = -0.05298011854;
		x[2*n] = a*x[deux_n_moins_un] + x[2*n] + a*x[2*n+1];

	}

	for(int n=0; n < p/2; n++){

		deux_n_plus_deux = 2*n+2;
		if(deux_n_plus_deux > p-1) deux_n_plus_deux -= p;

		//Prediction 2
		a = 0.8829110762;
		x[2*n+1] = a*x[2*n] + x[2*n+1] + a*x[deux_n_plus_deux];

	}

	for(int n=0; n < p/2; n++){

		deux_n_moins_un =  2*n-1;
		if(deux_n_moins_un < 0) deux_n_moins_un += p;

		//MaJ 2
		a = 0.4435068522;
		x[2*n] = a*x[deux_n_moins_un] + x[2*n] + a*x[2*n+1];

	}

	for(int n=0; n < p/2; n++){

		//Mise a l'echelle
		a = 1/1.149604398;
		x[2*n] /= a;
		x[2*n+1] *= a;

	}

	for(int n=0; n < p/2; n++){
		temp[n] = x[2*n];
		temp[n+p/2] = x[2*n+1];

	}

	memcpy(x, temp, p*sizeof(double));

	for(int i = 0; i< p; i++) {

		// std::cout << "Coef num "<<  i << " : " << x[i] << std::endl;
		if(file != NULL) fprintf (file, "%e\n",x[i]);
	}

	delete[] temp;
}


void synthese_97_lifting(double* x, int p, FILE* file){

	double a = 0;
	int deux_n_plus_un, deux_n_plus_deux, deux_n_moins_un;
	double * x_tri = new double[p];
	//double * y = new double[p];
	//y[0] = 0; y[1] = 2; y[2] = 4; y[3] = 6; y[4] = 1; y[5] = 3; y[6] = 5; y[7] = 7;

	//Mise en forme
	for(int i = 0; i< p/2; i++) {

		// std::cout << "Coef num "<<  i << " : " << x[i] << std::endl;
		x_tri[2*i] = x[i];
		x_tri[2*i + 1] = x[i+ p/2];
	}

	for(int n = 0; n < p/2; n++){

		//Mise a lechelle
		a = 1.149604398;
		x_tri[2*n] = x_tri[2*n] /a;
		x_tri[2*n + 1] = x_tri[2*n +1] *a;
	}

	for(int n = 0; n < p/2; n++){

		deux_n_moins_un =  2*n-1;
		if(deux_n_moins_un < 0) deux_n_moins_un += p;

		//Mise a jour 2
		a = -0.4435068522;
		x_tri[2*n] = a * x_tri[deux_n_moins_un] + x_tri[2*n] + a * x_tri[2*n+1];
	}

	for(int n = 0; n < p/2; n++){

		deux_n_plus_deux = 2*n+2;
		if(deux_n_plus_deux > p-1) deux_n_plus_deux -= p;

		//Prediction 2
		a = -0.8829110762;
		x_tri[2*n + 1] = a * x_tri[2 * n] + x_tri[2*n +1] + a * x_tri[deux_n_plus_deux];
	}

	for(int n = 0; n < p/2; n++){

		deux_n_moins_un =  2*n-1;
		if(deux_n_moins_un < 0) deux_n_moins_un += p;

		//Mise à jour 1
		a = 0.05298011854;
		x_tri[2*n] = a * x_tri[deux_n_moins_un] + x_tri[2*n] + a * x_tri[2*n+1];

	}

	for(int n = 0; n < p/2; n++){

		deux_n_plus_deux = 2*n+2;
		if(deux_n_plus_deux > p-1) deux_n_plus_deux -= p;

		//Prediction 1
		a = 1.586134342;
		x_tri[2*n + 1] = a * x_tri[2 * n] + x_tri[2*n +1] + a * x_tri[deux_n_plus_deux];


	}

	for(int i = 0; i< p; i++) {

		// std::cout << "Coef num "<<  i << " : " << x[i] << std::endl;
		if(file != NULL) fprintf (file, "%e\n",x_tri[i]);
	}


	memcpy(x, x_tri, p*sizeof(double));
	delete[] x_tri;



}
